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Ten Years in the Making — AUSM- family 


Meng-Sing Liou 

National Aeronautics and Space Administration 
Glenn Research Center 
2 1 000 Brookpark Road 
Cleveland, Ohio 44135 


We begin by describing the motivations that gave birth to the original AUSM scheme 
and then focus on the ingredients that has spurred its growth and acceptance by the wor 
of computational fluid dynamics. As it has played out more in the field, weaknesses have 
also surfaced. Hence, nutrients and supplements are prescribed to help it grow and stay 
strong and robust. In this paper, We will describe the saga of efforts owing to researchers 
who have contributed to building up the AUSM-family for the CFD community It is 
hoped that a healthy scheme will contribute to the accurate and robust solution of prob- 
lems encountered in a wide range of disciplines. We analyze numerical mass fluxes with 
an emphasis on their capability for accurately capturing shock and contact discontinu- 
ities. We will present a new scheme for the pressure flux, along with results for a host of 
test problems. 


Introduction 

C ONSIDERABLE progress in CFD has been made 
in solving equations of conservation laws over the 
last two decades, particularly in devising accurate and 
robust schemes for capturing shock and contact discon- 
tinuities. The ability to predict shock and cont act dis- 
continuities can be considered a prerequisite for a reli- 
able and accurate solution to both inviscid and viscous 
problems. The 1980s witnessed an explosive interest 
and research in upwind schemes for their capability of 
achieving high accuracy over a wide range of problems 
described by Euler or Navier-Stokes equations. Today, 
upwind schemes undoubtedly have become the main 
spatial discretization techniques adopted into nearly 
all major research and commercial codes. Yet some 
deficiencies or failures have been experienced by some 
upwind schemes, such as shock instability in multidi- 
mensions, creation of traveling waves in slowly moving 
shock, 1 and violation of positivity-preserving. 2 As 
CFD is being used more routinely and extended to 
more complicated systems of flow equations, the need 
for maximizing accuracy, efficiency and robustness for 
a wide variety of problems still remains the foremost 
concerns. Hence, the quest for the ultimate numerical 
flux scheme continues. 

Since the inception of the AUSM scheme in 1990, 3 
it has been adopted by researchers worldwide. It 
has been proven to be accurate, simple, robust , and 
easy to extend to other types of conservation laws, 
thus providing an attractive alternative to the existing 
schemes. In spite of the enormous progress achieved, 
deficiencies have been experienced, typically the post- 
shock overshoots and pressure oscillations along the 
transverse direction in the boundary layers, as sum- 
marized previously. 4 Several attempts have been made 
over the last ten years, e.g., Refs. [1,4-10] to improve the 


original scheme 3 and in general, some successes have 
been achieved. Now various versions of the AUSM-family 
schemes have been incorporated into both research and 
commercial codes. In this paper, we will show some of 
those results, including low Mach number flows, mul- 
tiphase flows, and DNS calculations. 

As evident in this paper and others, the numerical 
inviscid flux (which for the sake of simplicity we will 
hereafter refer to as numerical flux) plays a central role 
in effecting the success of a calculation, especially with 
regard to robustness and accuracy. Furthermore, we 
show in Ref. 2 that the mass flux plays the central role 
in the construct ion of a robust and accurate numerical 
flux that is simultaneously free of anomalies such as 
the odd-even decoupling and ” carbuncle’ 1 phenomena. 
This become cleat by realizing that the mass flux is 
common to the convective part of every" conservation 
equation of the fluid flows. 


This paper is organized as follows. First, we will 
take a broad approach to constructing the AUSM 
schemes, beyond the original one. We will examine 
in turn the mass flux and the pressure flux. Next, 
we will present the recently-introduced concept of a 
numerical speed of sound, which allows for a unified 
formula valid for the entire speed regime. It also lends 
itself convenient \y to the extension of the schemes to 
deal with multiphase/multifluid flows. Examples of 
applying the AUSM-family schemes to various types 
of flows will be shown. Finally 7 , we will propose a fur- 
ther development in regard to pressure flux, along with 
validation tests to demonstrate its effectiveness. 


N AS A/TM — 200 1-210977 


1 



Equations of Conservation Laws 

A set of equations of general conservation laws is 
considered: 


Qt k) ~ div(F , * ) - F^) = S. (1) 

We will denote by an overhead arrow “ the vectors 
associated with the Cartesian coordinates in three di- 
mensions. The conservative variables are given in Q (A) 
where the superscript *(= 1,2, ...) is introduced to in- 
clude multifluid models, often adopted for describing 
multiphase flows. 11 12 The inviscid and viscous fluxes 
are denoted respectively by F {k) and F v ( *\ whose 
definitions are omitted herein since they are rather 
standard. However, the source terms are dependent 
upon the physical problems studied. For multiphase 
flows , they can contain terms describing interfacial bal- 
ances of mass, momentum, and energy transfers due 
to phase differences /changes. We include this option 
because examples will be given later in the paper. 

During the 90s, a great deal of interest has been 
focused upon the development of a (local) precondi- 
tioning method to improve the convergence rate in the 
low Mach number regime. This is accomplished by 
premultiplying the time-derivative term with a condi- 
tioning matrix P 


rQ'*’ + div(F u ' ) - Fjt'J) = S, Q<*> 


{/;*.. f,.n.) r 




Several forms of the local preconditioning matrix T 
have been proposed in the literature, e.g., Refs. 13-1G. 


The discretization of viscous terms is rather stan- 
dard and is generally done with centered schemes. On 
the other hand, treatment of the source terms varies 


considerably and this subject is not so easy because 
it is quite problem-dependent and the terms can be 
extremely complicated. See for example the ones in- 
volved in the fluidized bed. 17 The subject is beyond 
the scope of this paper and will not be dealt with 
here. We shall restrict ourselves instead to the numer- 
ical representation of inviscid fluxes, which has been 
a subject of intensive effort by many researchers over 
the past three decades. 


Numerical Flux: AUSM-Family 

How It Began 

The 80s saw the rise of two classes of upwind 
schemes, namely the flux-vector and flux-difference 
schemes. The van Leer scheme belongs to the former 
and is perhaps the most robust one among all for the 
Euler equations, aside from a slight disadvantage in the 
shock resolution compared to the latter. Especially, it 
endures the positivity test and is free from shock in- 
stabilities, in addition to its algorithmic simplicity and 
generality. Interestingly, these properties served well 
during the era of revived interest in the hypersonic 
flight, e.g., the NASP program in the US and similar 


programs in the other countries. As more confidence 
has been gained by the CFD community in dealing 
with complexities in flow fields as well as in geome- 
tries, CFD has flourished and naturally Navier-Stokes 
solutions now have taken the center stage. The paper 
by A an Leer et al. lh points out that the flux- vector 
schemes are diffusive for the Navier-Stokes calcula- 
tions, as illustrated in the ID conical viscous flow 
where an incorrect wall temperature was predicted 
along with a thicker boundary layer. See Fig. (la). 
This spelt the downfall of the flux- vector schemes. 

In an attempt to resuscitate his scheme for t he 90s. 
\an Leer injected the flavor of flux-difference split- 
ting, as first suggested by Hanoi and Scbwane, 19 in 
the CFD symposium held at. NASA Lewis Center in 
1990. 1 The improvement due to this new idea is clearly 
demonstrated in Fig. (lb), showing the wall temper- 
ature close to the correct value of 13.7 ( Pr ~ 1.0), 
but unfortunately afflicted with pressure irregularity 
at the edge of the boundary layer. It appeared that 
the flux-vector scheme was a phoenix, the catalyst for 
its rebirth being the combination with flux-splitting 
scheme. So Van Leer asked the question: Can pure 
flux- vector splitting be saved? The key word is ‘’pure” 
and the answer may be still up in t he air. But the ques- 
tion can be rephrased as: Can the Van Leer splitting be 
engineered so that good “genes” are kept? The answer 
is very likely. Hence, research started and the original 
AUSM scheme 3 (presented in this conference the first 
time in 1991 20 ) began to take shape. As a result of 
this quest, the AUSM scheme gives the temperature 
Profile in excellent agreement with the solution by the 
Godunov scheme. 21 


How it is shaping up 

In this section, we will look at the algorithm in- 
volved in the AUSM-family schemes and the new de- 
velopments. For more details and other numerical 
properties, the reader should consult, with the cited 
references. 

As a first step common in the AUSM schemes, we 
explicitly split the inviscid flux (written in three di- 
mensions) into two parts: 

F = Y {c) |P = m^-bP, (3) 

where 

m = p\ , jj; = (U u, t\ tc, H) t , P = (0, pi,pj,pk y 0) T . 
(4) 

: This occasion brought together the fathers of two flux- vector 
splittings, Joe Steger and Bram van Leer, although Joe talked 
about an entirely different subject--the chimera method. My 
family and I had a great deal of fun with them in my home that 
evening, the fun of course was heightened with Bram’s playing 
piano, especially performing Joe’s favorite piano and orchestral 
piece, Symphonic Variations by Cesar Franck. 
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Van Leer F.V.S. 


Van Leer/Hanel91 




a) Van Leer FVS. b) Van Lcer/Hanel (’90). 


i-Splitting) 



c) AUSM (’91). 


Exact Riemann Solver 



d) Godunov Riemann 
solver. 


Fig. 1 Hypersonic conic flow, M^c — /.95, 0 rO ne/2 — 
10° and Reoo = 4.2 x 10 5 . 

The first term in F is the convective flux F (c) . indi- 
cating the convection of ip by the juass flux tn and the 
second term is the pressure flux P. containing nothing 
but the pressure. It is noted that the conservation of 
the total enthalpy is guaranteed if H , instead of to- 
tal internal energy (£), is contained in ip. Halt and 
Agarwal 22 split H into E and pj p and put them re- 
spectively in the ip and P terms and called it the WPS 
scheme. 

In terms of the component fluxes in the directions 
of x, y, and z, we have 

F = (ro,i!i +P x )i"t - P y )J -r (m-t/> + P z)k. (5) 


where 

ni[ — pui , P / = (0, p&ix, P&izi 0) 5 l = 

( 6 ) 

Tn a control volume, the mass flux m n through a 
control surface element having a unit normal vector 


fi = (n x ,n y ,n z ) is given by 

m n = pV ■ n = n x m x ~rTi y m y —n~ ni~ . (7) 

And the associated flux becomes, 

F„ =m„ i' -r P», P„ = p(0, n x , n v , n : ,0) T . (8) 

Formally, this equat ion looks the same as that along 
an individual Cartesian coordinate direction. Hence, 
at each control surface, the mass flux is treated in an 
one-dimensional fashion. At the discrete level, this 
is also what one needs to do for defining the flux at 
the interface in a finite volume. Hereafter, we will as- 
sume that this local orientation has been accomplished 
and the velocity vector (hence mass flux) has been de- 
composed into components normal and parallel to the 
surface vectors fi. Therefore, the subscript n denot- 
ing the normal component will be dropped. 

From the above equations, the principal quantities 
in F u are again the convective flux and the pressure 
flux. The distinction of these two fluxes gives rise 
to the basis for the development of the AUSM-family 
schemes. Since the mass flux appears in all equations, 
its effects will be felt by all the variables. Hence, we 
believe that it is desirable to observe this fact at. the 
discrete level as well when devising a new scheme. Sig- 
nificant. benefits can be derived as well. For example, 
the numerical dissipation term is scalar even for the 
system of equations; it is just as easy to add more 
conservation equations insofar as the numerical flux is 
concerned. 

It is possible to writ e a numerical flux, mimicking 
the expression at the continuum level, in terms of a 
common mass flux in the following general upwind 
form. 

ft/2 =Wl/2 1pL/R + Pl/2 

= 77 ly/ 2 VL - m \/2 — Pl/2- V ’ 


Here the contributions of and i!)r are weighted by 
the split masses (?n 1/2 ,m 1/2 ), which must follow the 
consistency requirement, 

m,/ 2 = m + l/2 -f- m 1/2 - (10) 

This fact is automatically sat isfied by the first, element 
of f. One can rewrite Eq. (9), using Eq. (10), as 

ft/2 = ^ mi/ 2 (i’Ld-if’«)-^X , m (^fl-^i.)+Pi/2i (11) 

The dissipation term, V m , is 

V m = m l/2 -m l/2 . (12) 

The subscripts “L” and “R” are understood to mean 
the cell centers on either side of the interface at which 
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the normal vector is assumed to point from “L” to 

“R” . 

The quantities (wii/ 2 > Wi i/ 2 ) are required to satisfy 
these conditions, 

• + • - 

( m iA>) > 0, (w»i /2 ) < 0, (13) 

so that they provide proper upwinding, thus ensur- 
ing stability. In the AUSM-familv schemes, these two 
variables are mutually exclusive, i.e., 


For the AUSM + : we first define the interface Mach 
number 

(18) 

then 

\ /2 - |A/,/ 2 |). (19) 

Now for the AUSMDV: we first define the interface 
split Mach numbers, 


(m l/2 ) (m 1/L> ) =0, (14) 

It must be noted that the flux expressed in the form 
°f Eq. (9) implies that the numerical dissipation is 
of the scalar, rather than the matrix form, because 

the same factors m x /2 are applied throughout for all 
conservation equations. The flux difference splitting 
schemes are known to belong to the category of matrix 
dissipation. On the other hand, the category of scalar 
dissipation encompasses several existing schemes other 
than the AUSM-fainily schemes, such as central dif- 
ferencing with artificial damping, the Van Leer flux 
vector scheme, the HLLE scheme. Indeed, there are 
several attractive properties associated with the scalar 
form of dissipation. From the algorithmic viewpoint, 
it offers simplicity, efficiency, and generality allowing 
for an easy extension to other systems of equations. 

The 1991 AUSM scheme has served well by laying 
out the basis for further developments. One of the im- 
portant developments is the concept of common speed 
of sound, which makes an accurate resolution of con- 
tact and shock discontinuities possible for both steady 
and unsteady flows. As a result, two new members 
of the AUSM-familv were generated, We shall in this 
paper, specifically concentrate on a unified formula- 
tion encompassing both the AUSM+ and AUSMDV 
schemes. 

In what follows, we will give some basic formulas 
used to define the mass flux. To facilitate the discus- 
sion, we first define the following split functions. 


together with the blending functions, 


■+ 


-7, 


. / = p/p- 


~ 1/2 /,-/« 

Then, the interface Mach number is 


( 20 ) 

( 21 ) 


M i/ 2 = A/+ 2 + Jlff /2 . (22) 

It is noted that the construction of the interface split 
Mach numbers M+ 2 in the AUSMDV scheme is some- 
what similar to that in the Van Leer’s flux scheme, 
one might then wonder if the AUSMDV would be also 
afflicted with the same shortcomings. To the cred- 
its of the variables. varying with flow variables, 
they make AUSMDV an accurate scheme for capturing 
contact discontinuities, hence appropriate for viscous 
solutions. 

We stress again that a common speed of sound 
a i /2 — U*) is used in the formulation in defin- 

ing the U L and l Mach numbers and in Eq. (9). 
In Ref. 6, we give a special formula for a x/2 so that 
an exact capturing of a stationary normal shock is 
achieved. Ot herwise, any averages of the “L” and “R” 
states should be appropriate. More importantly, this 
possibility of flexibly allowing other definitions of the 
common speed of sound opens a very rewarding op- 
portunity, as will be discussed later. 


£>< 


M%.(M) = ±UM± l) 2 


(15) 


Vow the mass flux is immediately available by using 
the quantities 




if \M\ > 1, 
( 1 — 2.Vf 7), otherwise, 


m 


- a\rApi. A /2 -j- PhM 1/2 ) 

= /-iWx/Apl - Pr) - D p (p R ~ p L )] 


and 


nUfA/) = 




>^[(±2 - M) T 3A/.A4*,], otherwise 


(IG) remark that a clear difference between the AUS- 

MDV and AUSM+ schemes, insofar as rn is concerned, 
is in the definition of 

if | A/ 1 > 1. eaS y i 0 s h ow that: 


(17) 

The numerals in the subscript of 

D - J 

I AUSM+ 

- M~ 2 AUSMDV 

(24) 

and yy indicate the degree of polynomials. ^ 

Using a common speed of sound a 1/2 to define M L = 
u 'u/ a i /2 and Mr = u„ B /a 1/2 . 


= 0 AUSM+ 

< 0 AUSMDV 

(25) 
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The dissipation coefficient in the convective flux, in 
Eq. (11), is now explicitly given as 

D m = I m l/2 I 

We can also rewrite the mass flux in the most general 
form, 

m,/ 2 = (m) - -^X>(Uz,,Uft). (27) 

where (m) is a sort of centrally weighted aver- 
age, but the detail is not important for our discus- 
sion. However, we shall focus on the dissipative term 
X>(U l,U r ) which can be further expanded in terms 
of differences of primitive variables U - {p. V.p) as 
follows. 

V = V {p) {TS)&p + ^D (u,) (U)Au, - P (p) (U)A/>, 

(28) 

where U(U/.,U fl ) are some mean quantities, and 
the difference operator is A(*) = (•)/? - (•Ji- 
lt. is shown in Ref. 2 that the fact whether the 
pressure dissipation coefficient D^ vanishes for all 
conditions plays a decisive role in determining the oc- 
currence of the anomalies in shock instabilities. 

Finally, another important variable is the pressure 
flux, which may be written as 

P„2 =P 1/2 (0,ffi £ ,n i/ ,n J ,0) 7 . (29) 

Clearly, all one needs in the pressure flux p 1/2 is simply 
the definition of Pi/o* 

In all the AUSM-family schemes, the interface pres- 
sure has been simply given by 

p 1/s = V+ } (Ml)P:. + ^(5, ( M «)p* • ( 3 °) 

As simple as it may seem, there apparently are 
enough opportunities to further enhance the AUSM- 
family. A new version of the pressure flux will be 
presented later. 

The accuracy of the AUSM"*” scheme was thor- 
oughly established by Darracq et al. in Ref. 23 
in which they have done studies of grid refinement 
and spatial order of accuracy for several airfoil flows, 
against the measured data. Tables 1 and 2 present 
the comparisons of calculated results for 2D and 3D 
turbulent flows. They concluded that “for all runs the 
AUSM+ predictions agree better with the experimen- 
tal data than results obtained with the Roe scheme.” 
In general, the differences between predictions from 
these two schemes become more apparent in the pre- 
diction of flows near the leading edge on the suction 
side. It also shows that AUSM+ solutions converge to 
the grid independent solution faster than those of the 
Roe splitting. 

These results, along with those from my own and 
others, seem to suggest that the ALSM schemes (es- 
pecially AUSM+ and AUSMDV) yield little numerical 


Table 1 Comparison of lift and drag coefficients for 
the RAE-2822 airfoil, AUc =0.73, Q— 2.79°, Rene = 
6.5 x 10 6 , Baldwin- Lomax model. (Ref. 23) 


Scheme 

Order 

Mesh 

C,. 

Cd p 

Cd tot 

Roe 

2nd 

coarse 

0.7755 

0.0153 

0.0209 

AUSM+ 

2nd 

coarse 

0.7931 

0.0144 

0.0200 

Roe 

3rd 

coarse 

0.7814 

0.0136 

0.0192 

AUSM+ 

3rd 

coarse 

0.79G1 

0.0133 

0.0188 

Roe 

2nd 

fine 

0.7916 

0.0133 

0.0189 

AUSM+ 

2nd 

fine 

0.8046 

0.0132 

0.0187 

Roe 

3rd 

fine 

0.7927 

0.0131 

0.0187 

AUSM+ 

3rd 

fine 

0.8064 

0.0131 

0.0187 

Expt." 



0.803 


0.0168 

Table 2 C 

Comparison of lift and drag coefficients for 

the ONERA-M6 airfoil, M 

^.=0.84, 

o=3.06 

\ He.oo = 

1.749 x 10 7 

Baldwin-Lomax model. 

(Ref. 23) 

Scheme 

Order 

Mesh 

Ci. 

C Dv 

C'Dtot 1 

Roe 

2nd 

coarse 

0.2558 

0.0168 

0.0220 

AUSM+ 

2nd 

coarse 

0.2655 

0.0164 

0.0215 

Roe 

3rd 

coarse 

0.2604 

0.0141 

0.0192 

AUSM+ 

3rd 

coarse 

0.2667 

0.0138 

0.0189 

Roe 

2nd 

fine 

0.2782 

0.0137 

0.0188 

AUSM + 

2nd 

fine 

0.2819 

0.0138 

0.0188 

Roe 

3rd 

fine 

0.2791 

0.0132 

0.0183 

AUSM+ 

3rd 

fine 

0.2825 

0.0132 

0.0182 


dissipation. Hence attempts have been made recently 
for simulations touted as demanding high accuracy, 
such as Large Eddy Simulation (LES), and Diiect Nu- 
merical Simulation (DNS). 

Recently, Billet and Louedin 25 combined the AUSM 
scheme with a very interesting adaptive limiter (named 
triad limiter, triad) to gain high accuracy for DNS- 
type simulations of unsteady flows. The limiter is still 
built upon the third-order accurate MUSCL formula- 
tion, but the accuracy of the results rivals that of the 
higher-order schemes. Advect.ion of a Taylor vortex is 
simulated using the AUSM-^.w scheme 25 and the 
velocity and pressure profiles, shown in Fig. 2, are 
in excellent agreement with those from a sixth-order 
accurate Hermitian scheme. Also Fig. 3 displays the 
time development of a 3D mixing layer using the same 
scheme, show’ing a nearly identical result as that from 
the DNS. 

Other DNS/LES calculations, for example, can be 
found in Ref. 26-28. There appears a common strategy 
in the simulation of low speed flows, insofar as using 
the AUSM schemes is concerned. That is, the pressure 
term is modified by replacing it with a simple average 
of neighboring pressures. This seems to have worked 
well, giving smooth solutions without even adding nu- 
merical dissipations. However, this simple replacement 
gives rise to a smeared shock when applied to super- 
sonic flows. Nevertheless, what this indicates is that 
there are still some things to be done in the area of 
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f. ,g ^ 2 - Ta ^!° r ™ rtex advect ion: longitudinal distribution along the centerline of variables 

(bottom). The blow-up view near the extrema are shown on the right. (Ref. 25) 


v (top) and p 


pressure splitting. This is the new development yet to 
be disclosed later in this paper. 

In what follows, we will discuss another important 
chapter about the AUSM family. This is the extension 
to the low Mach number flows and something other 
than the aerodynamic flows, for example, multiphase 
flow's. 

As it has become known during last decade that 
the detrimental deficiencies in forcing the compress- 
ible upwind codes onto solving low' speed flow's are : 
(1) extremely slow or stalled convergence, more so as 
the flow speed decreases and (2) the flow solutions 
can be globally incorrect (rat her than just locally as 
in the case of smearing shocks). These two phenom- 
ena are not related because the first one originates at 
the continuum level, depending on the form of govern- 
ing equations being solved, irrespective of whether the 
scheme is centered or upwind. However, the second 
one is inherently tied to the upwind scheme where the 
eigenvalues, strongly depending on the usage of the 
speed of sound, are employed. 

In the 90s, active research has been conducted to 
conquer the first problem, in the name of local pre- 
conditioning, such as those by Van Leer et ah, 14 
Turkel, 13 and Merkle, 15 and their subsequent publica- 
tions. While there are differences in approaches, they 
all attempt to achieve the same objectives of making 
the eigenvalues of the new system of equations, Eq. 
(2), the same order of magnitude. A condition num- 


ber k defined as the largest ratio of eigenvalues, 


■ oc, as fa 


■ 0, and a held fixed. 


is a useful measurement.. Clearly there is a large dis- 
parity of wave speeds as |u| -» 0 and as a result, this 
has been identified as the source of slow (or no) con- 
vergence. 

In the preconditioning strategy, one can think of 
seeking to modify the system in such a way that the 
corresponding speed of sound would be altered to be- 
have like |u| as it approaches zero. 

Consequently, we will define the numerical speed of 
sound by 

d = /(A/;A/»)a, (32) 

where the scaling factor may be of this form, 


/(M; A/,) = 


(1 - M?)h\P ~ 4M'i 
1 + A/? 


and the reference Mach number, 

= min(l, max(A/ 2 , A/ 2 0 )). (34) 

The cutoff parameter M co is introduced to prevent a 
singularity at stagnation point. It is a user-specified 
parameter and A/ 2 0 = Hr 4 has been used. Details can 
be found in Ref. 29. 

Now the condition number becomes, 

|tt| -f- CL 

“ , | ~ -+ 0{ 1), while a -> 0, as |u| -*■ 0. (35) 
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a) Time evolution of vorticity thickness. 



b) One static pressure surface at t=45. 

Fig. 3 3D mixing layer. (Ref. 25) 

That is, the condition number remains order of unity 
at low speeds. The numerical dissipation based on this 
new speed of sound now scales with the local speed |u|, 
instead of the local speed of sound a. As a result, the 
accuracy can be restored as it is applied to low Mach 
flows. 

Now we can define the new Mach numbers based on 
this numerical speed of sound a as 

Ml/r = ^ (36) 

and these would be the entries to the equations for 
A1(4) and 'P ^ . This version is then denoted with suffix 
a. such as AUSM + -a. 

Now we show the results performed by Mary et 
al. 30 using the low Mach number version of the 



a) Snapshot of vor- 
ticity. 



2 


o-i— . ' 

0 25 50 75 100 

time 

b) Time history of vorticity 
thickness. 


Fig. 4 Comparison of solutions for subsonic com- 
pressible mixing layer, third-order accurate so- 
lutions of three different time step sizes (At = 
0.03,0.06,0.4 (shown by ...) are compared on the 
same number of grid points with that by the sixth- 
order Hermit ian scheme (denoted by □□). Results 
from two small time steps are indistinguishable. 
(Ref. 301 



a) Initial setup. 



•0 ' 


*1 

b) After interaction. 


Fig. 5 Interaction of a Gaussian temperature spot 
with a shock. Comparison with the 4th order accu- 
rate WENO scheme (denoted by □□) with the 3rd 
order AUSM+ solutions with three different values 
of the limiter compression parameter ©. (Ref. 30) 


AUSM + scheme with a third-order accurate spatial 
interpolation. Figure. 4 shows that the results of 
a subsonic compressible mixing layer are in excellent 
agreement with a sixth-order Hermitian scheme on the 
same number of grid points, hence demonstrating the 
accuracy of the AUSM+ scheme. Another example, 
given in Fig. 5, compares the result with that of the 
4th-order accurate WENO scheme for a temperature 
spot interacting with a shock, again revealing the high 
accuracy of the AUSM+ scheme. 

Implementing AUSM + -a in the code results in a 
significant improvement not only in solution accu- 
racy, as seen above, but also in convergence rate as 
well. Another unforeseen but pleasant consequence 
is that the pressure oscillations, observed along the 
transverse direction in the viscous layer when using 
AUSM + scheme, are no longer there. This results 
from the fact the much reduced numerical dissipation 
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Fig. 6 Pressure contours for the shuttle external 
tank problem for M = 0.01. a): using the standard 
AUSM+ at N=6400 time steps; b): using numer- 
ical speed of sound at N=1000 time steps, and c): 
magnified view near the nose. 

now scales properly with the pressure variations, as 
seen in Fig. 6. 

Another case in point is the application to a 3D 
low speed flow over a high-lift three-element trap- 
wing configuration. The flow conditions are: M <*, = 
0.1498, = 14.7 x 10 6 and angle of attack of 20 

degrees. The computational geometry model con- 
sists of a body pod, a wing, a full-span slat, a full- 
span flap and the tunnel walls, as displayed in Fig. 
7. Rogers et al. performed an extensive numerical 
study of the aerodynamic characteristics of this con- 
figuration using the preconditioned version of the Roe 


Fig. 7 Trap wing model in a wind tunnel, the 
tunnel grid is plotted every fourth grid point. 

scheme in the OVERFLOW 31 code. The effectiveness 
of the AUSM+-a scheme was tested for this configu- 
ration and details will be given in a separate paper. 32 
The calculation 2 was also performed with the OVER- 
FLOW code, with the Spalart-Allmaras one-equation 
model. 33 The pressure distributions at various span- 
wise locations are compared with the experimental 
data in Fig. 8 and they are in excellent agreement. 

The vortices generated from the wing tips and the 
body pod are illustrated with the particle traces dis- 
played in Fig. 9 

Recently, the AUSM-family has been extended to 
the multiphase flow calculations, e.g., in Refs. 34-36. 
Paillere et al. solved a system of two-fluid models 
with interfacial source terms included. Several fea- 
tures that are different from the usual equations for 
aerodynamic flows add complexity significantly. That 
is, the system is no longer in conservative form because 
of the presence of the source terms and the system 
is not guaranteed to be hyperbolic because it admits 
complex eigenvalues. Figure 10 displays the computed 
evolution of an initially homogeneous mixture of liq- 
uid water and air, under the action of gravity, moving 
toward a complete phase separation at the final steady 
state. The phase separation begins at the both ends 
and gradually migrates towards the center, as illus- 
trated by the evolution of void fraction of the air, with 
x measured from the top; the pressure at the steady 
state is essentially constant in the gas region, but in- 
creases linearly in the liquid region, as it should. Next, 
Paillere et al. also obtained results for a water/air 
column oscillating under the gravity. The solution in- 
cludes the effect of interfacial drag. The void fraction 
of air and the time variation of liquid velocity at the 

2 L)r. Stuart Rogers of NASA Ames Research Center, kindly 
provided the grid, the post- processing code for extracting the 
surface pressure and consultations. 
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1.0 2.0 3.0 

X/C 


Fig. 8 Pressure distributions. 

bottom of the tube are shown in Fig. 11. displaying a 
cyclic motion. The initially sharp profile of void frac- 
tion is now smeared at t = 20, due to the interfacial 
drag. 

Another example of multiphase flows involves a wa- 
ter flow over a hemispherical cylinder. The flow can 
undergo cavitation if the pressure difference ( cavita- 
tion number, K = 2(poo - p v )/p ocU^) is low enough. 
Figure 12 shows the water density contours in a tur- 
bulent flow under various cavitation conditions. It 
clearly shows the phase transition between the liquid 
and vapor states. As K decreases, the pressure in the 
expansion region drops to the vapor pressure, resulting 
in the generation of a vapor phase and the growth of 
a cavitation bubble. Pressure recovery further down- 
stream leads to the collapse of the cavity in a “wake” 
region. The structure of the wake region is strongly 
influenced by both the thermodynamic model and the 
velocity field, which in itself is influenced by the tur- 
bulence model. 

Another recent accomplishment of using the AUSM 
scheme has been reported in Ref. 17 by De Wilde et 
al. for an extremely complicated set of equations and 
flow patterns involving solid particles and gas phase in 
an industrial scale riser with a diameter of 1 .56 m and 
a height of 14.434 m. The system is described by an 



a) Particle traces from the body pod. 



b) Particle traces from the wing. 


Fig. 9 Particle traces. 

unsteady 3D turbulent, two-fluid model. The source 
terms include gravity, buoyancy, and stresses due to 
gas-solid interactions. The solid particles along with 
the gas, entered at velocities of 6.0 and 12.635 (m/s) 
respectively, in an inlet at the bottom of the riser, and 
the mixture exited at the top. Due to the inelastic 
particle-particle collisions, flow instability is triggered 
and a periodic slugging flow pattern was obtained by 
them. Figure 14 displays the evolution of solid volume 
fraction in one cycle. A perturbation is seen to origi- 
nate at the top of the riser, grow underneath towards 
the bottom of the riser, and reach the maximum extent 
at about 3.2 sec and then move upward till it is blown 
out at the top when t=6.4 sec. Then the cycle contin- 
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a) Physical description. 
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b) Evolution of void c) Pressure profile at 
fraction to steady state. steady state. 

Fig. 10 Phase separation test case. (Ref. 34) 


ues. From the results, one sees a large scale motion in 
the axial direction together with a radial variation in 
each time frame. However, no circumferential asym- 
metry was found for the conditions calculated. The 
authors state that the oscillation frequency of 0.15 Hz 
is in good agreement with that reported in literature. 
The time-averaged result give distinct boundaries of a 
cell-like structure. 

The above has merely captured representative ad- 
ventures of the AUSM-familv schemes beyond aero- 
dynamics. The DNS/LES calculations demonstrate 
the accuracy inherent in the AUSM-family schemes, 
rivaling that of the higher-order schemes. The multi- 
phase flow problems are certainly far more difficult 
to deal with, not only from the closure (modeling) 
point of view, but also from the algorithmic one. The 
source terms strongly couple variables associated with 
all phases and give rise to an extremely stiff system. 
Robustness of a numerical scheme is the key to the 
capability of simulating these types of flows. 

To add a contribution for the new millennium, I 
shall present in wdiat follows the result of a recent ef- 
fort, with an aim at further improving the scheme’s 
robustness and accuracy. Up to this point, calculating 
flows at low speeds often required adding a pressure- 
diffusion term to the mass flux in order to enhance 
convergence. This has been commonly done in the in- 
compressible code to ensure pressure- velocity coupling. 
In spite of its effectiveness, it nevertheless stems purely 
from numerical consideration. So the questions to ask 
are: (1) whether this (adding the pressure- velocity 
term) is absolutely necessary? (2) if yes, whether this 
is the only way? and (3) if not, what then? 



Oscillating Manomalar 



a) Physical setup. 


b) Initial conditions. 


Oscillating Manometer 



Oscillating Manometei 



c) Void fraction at t— 20 d) Liquid velocity at 

s. x=10 m. 

Fig. 11 Oscillating manometer. (Ref. 34) 

p = 1061 kg I’m* P = 29 kg^m* 



ctntcrline 



K-O.fl (no ctvilaiion) 


Fig. 12 Density contours of liquid water flow over 
hemisphere/cylinder for various cavitation num- 
bers. (Ref. 35,36) 

The answer to (1), based on my experiences and 
those reported in the literature, is that it is desirable 
to have this sort of mechanism, although it may not 
be absolutely necessary. Thus, question (2) leads to 
finding an alternative, and the more difficult question 
(3) will be left alone for now'. 

In AUSMDV, 1 the convective part of the momentum 
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Fig. 14 Time-accurate calculations of the solid volume fractions in an industrial scale riser, with gas and 
granular particles flowing from the bottom and exiting at the top. (Ref. 17) 
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Fig. 13 Effect of spatial discretization accuracy on 
surface pressure distribution. (Ref. 35,36) 

flux consists of blending the flux- difference and flux- 
vector procedures (thus denoted with DV). A notable 
advantage of this strategy is that it gives smooth shock 
profiles, e.g., in shock-shock interactions. The notion 
was adopted in the AUSM+-W, 4 but its effectiveness 
has never been extensively tested. This blended flux 
can be recast and the extra terms can be reassigned to 
the pressure flux, giving rise to a new pressure flux con- 
taining a term proportional to the velocity difference. 
This interpretation, however, has a sound connection 
to the characteristic equations, 

dx 

dp ± padu — 0, along — — u ±a. (37) 

An integrated form for the interface pressure for |u| < 
a is 

Pl/2 = ^[(Pl ~Pr) ~ Pl/2 a l/2(u R -lij]. (38) 

Hence, we can beef up the pressure flux, Eq. (30) by 
including the velocity difference term, 

PI/2 = V+^M^Pl + V~ 5) {M R )p R 

-V+)( M L)V- 5) {M R ) Pll2 a\ /2 {M R - M l ), 

(39) 

where for the interface quantity pip 2 we may use, e.g., 

The coefficient involving and is 

introduced to automatically transition between super- 
sonic and subsonic conditions. The pressure now is 
explicitly coupled with the velocity field by the addi- 
tion of the velocity difference term. As M tends to 
zero, the pressure flux reduces to a form similar to Eq. 
(38), but with half of its coefficient, 

Pl/2 = \[{P,. +Pr) ~ ^Pl/2 a l/2( u n - «,.)]. (41) 


AUSN/r-u 



Fig. 15 Hypersonic conic flow. 


To denote this new version, we use the suffix U u” (for 
velocity diffusion) and call it AUSM+-U. Thus, when 
Tn. p 29 is also included, the scheme reads AUSM + -up. 
Or if the numerical speed of sound is also act ivated, the 
version becomes AUSM+-au. Note that if AUSM + -au 
is used, then the coefficient in Eq. (39) is scaled how- 
ever, with the numerical speed of sound a. This is 
just what we want for low Mach number flows since 
a = 0(u) as \u\ — > 0 and hence the coefficient is 
scaled by the magnitude of local velocity. More de- 
tails concerning this latest development shall be given 
in a separate paper. 37 

In what follows we shall consider several benchmark 
problems I usually used for testing numerical schemes. 
They represent various facets encountered in typical 
flow problems. We shall first consider one-dimensional 
problems and use the first-order scheme. 

First, we must require that the new pressure flux be 
capable of correctly predicting viscous flows, such as 
the hypersonic conic flow mentioned in the beginning. 
The new scheme, AUSM+-U, as in the other AUSM- 
familv schemes, gives the correct solution, as seen in 
Fig. 15. The reason that the velocity difference term 
does not cause adverse effect is that the velocity com- 
ponents in Eq. (39) are those normal to t he cell face 
and they are continuous across the viscous layer. 

The complaint voiced commonly about the 
AUSM+ concerns the overshoots resulting from 
strong shock-shock interactions. Figure 16 shows the 
comparison of results from various schemes. The new 
pressure flux AUSM+-U now gives a significantly im- 
proved result and the AUSM + -up scheme completely 
removes the overshoots, yielding results as good as 
those by the AUSMDV, Roc and Godunov schemes. 
It is a major success since the AUSMDV scheme in 
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c) Roe’s scheme. 


d) Godunov’s scheme. 


e) Roe’s scheme. f) Godunov’s scheme. 

Fig. 16 Colliding shock problem 


this regard. This test clearly indicates that the new 
pressure flux formula is a worthy replacement of the 
old one and will be used again in the following tests. 

The second problem concerns a shock moving slowly 
against a flow, as studied in Ref. 38. Figures 17 
show's the strength of linear and nonlinear waves by 
the AUSM+, AUSM+-U (AUSM+-up is indistinguish- 
able and not shown), Roe, and Godunov schemes. It is 
known that the Roe and Godunov schemes produces a 
noticeable long wave trailing the shock, as seen in the 
figure. The AUSM schemes however, perform quite 
well. 

The third problem is a shock moving through a con- 
stant area channel in which the grid at the centerline is 
perturbed alternately at odd and even points, as pro- 
posed by Quirk. 39 Figure 18 displays the result from 
the new scheme, it like the AUSM+ is clearly free of 
any shock instabilities. 

Next we will consider two two-dimensional super- 
sonic problems using a third-order spatial discretiza- 


Fig. 17 Slowly moving shock problem. 
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Fig. 18 Odd-even problem. 
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a) Pressure contours:AUSM + -au 




b) Residual history 


Fig. 19 Supersonic flow in a ramp-channel, M-x = 



10 ° 
to*-' 

10-2 
10" 3 
I 

10" 5 
10-* 

10" 7 
10-® 

0 2000 4000 6000 

N 


AUSM*-ou 



a) AUSM+-ap. b) AUSM+-au. 

Fig. 20 Residual history for the shuttle external 
tank at various Mach numbers. 


tion. A supersonic flow in a ramp-channel is shown 
in Fig. 19. The pressure contours show a smooth be- 
havior across shocks and the residual converges mono- 
tonicallv, again reaffirming the effectiveness of the new 
pressure flux. 

Figures 20 display the comparison of the conver- 
gence rates for different flow speeds between the two 
schemes. AUSM+-ap 29 and the present AUSM + -au. 
The convergence history by the new pressure flux is 
essentially the same as the AUSM + -ap scheme up to 
N = 2400, but appears to have slowed down afterward. 
The reason is not clear and will be further investigated 
using different codes and for different problems. How- 
ever, it must be noted that the solutions from these 
two schemes are not distinguishable. 

The final test is a blunt body problem often used 
by Radespiel. 40 This problem has several features to 



a) New pressure flux, b) Blow-up view. 

AUSM+-au. 

Fig. 21 Blunt body problem, M 0 c = 10. 


study. The grid 3 is clearly for the viscous calculations, 
but is tested here for the Euler solutions. First, the so- 
lution is free of carbuncle phenomena, consistent with 
the conjecture given in Ref. 2 since t here is no explicit 
pressure diffusion term in the mass flux. Secondly, the 
pressure contours are smooth, not. only near the wall, 
shown in the blow-up view near the stagnation region, 
but. also near the sonic line region. 


Concluding Remarks 

We have just seen a brief history of the AUSM- 
family, having finished the first ten years from its 
inception. It has seen both triumphs and setbacks 
while facing realities and difficult tasks. Nevertheless, 
many researchers have contributed in different ways 
towards its growth. The question is whether it will 
have the longevity to beat the odds in the future, to 
see a even wilder world. I believe that it has the fun- 
damentally right stuff, even though some turns and 
twists are expected every time there is a new fron- 
tier to be explored. In this brief tour, we have seen 
the successes, not only in aerodynamics, but also in 
the areas of DNS/LES and multiphase flows. Further- 
more, the new pressure flux boosts the strengths of the 
AUSM-family, having removed the overshoots behind 
the shock-shock interactions-a major success since the 
AUSMDY scheme in this regard. These should give us 
sufficient confidence in its ability to provide needed 
accuracy, efficiency, and robustness. 


3 The grid was provided by Prof. Rolf Radespiel, Braun- 
schweig Technical University, Germany. 
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